*------------------------------------------------------------------------------
*
*By:  David E. Lewis
*Date:  May 25, 2001
*
*My hope is to read the file in to stata, take care of missing data
*set the data for survival analysis, and estimate all of the models for my
*insulation and durability paper.
*------------------------------------------------------------------------------
*First, I need to do some housekeeping.
set memory 128000
set l 72
set level 95
set matsize 150
version 7

*Now, I'll look at the spell file.
use c:\stata\durdata51

*I am going to drop to courts in the data set, the Temporary Emergency Court of Appeals and the 
*Appeals Court for the Federal Circuit.  The only court-type agencies in there are in cabinet 
*departments, one for the armed services and the other for veterans affairs.

drop if agencyid==10439
drop if agencyid==10329


*Now I'll describe the data.
describe

*Now I will take care of some missing cases.
keep if cumdur~=.
keep if year~=.
keep if estate~=.

*cumdur is cumulative duration or age.  year is what is appears to be.
*estate is end state.  It is 0,1 depending on whether an agency is active
*or terminated.

*summarize
 
*First, I need to take care of one problem which is that if an 
*agency was terminated on the first of the year its cumdur is 0 instead of
*1.  I also need to take care of the problem that if an agency was
*terminated on the first of the year that its cumdur=cumdur for the 
*previous case.
replace cumdur=1 if cumdur==0
replace cumdur=cumdur+1 if cumdur[_n-1]==cumdur

*Now I will try to get stata to recognize the spell file.  stset
*tells stata this is a spell file, cumdur is the agency age in days
*id is the unique agency id code.  There are multiple observations on single
*agencies.  Scale tells stata to report the results in years.

stset cumyear estate, id(agencyid) 

*Now I will describe the survival-time data.

stdes
more

*Now I will summarize the data

stsum
more

*First I graph the product limit estimate of the hazard rate.  This isn't so helpful and 
*looks decreasing but non-monotonic.

ltable cumyear estate, hazard graph s(i) c(m) l1(Hazard Rate) t1(Appendix C. Product-Limit Estimates of Hazard Rate)
more

*I have to create some of the variables that will be useful in the analyses.
*I have to make the difference measure an absolute value.
generate abdistnc=abs(cshmed-cspres)
generate abhsecha=abs(hsechang)
generate abpresch=abs(preschan)

*Now I will generate the interaction terms.
generate partintr=unfrmaj*unfrpres*unified
generate ideointr=abhsecha*abpresch*abdistnc

*Now I will change the party measures to republican indicators
generate reppres=1 if dempres==0
replace reppres=0 if dempres==1
generate rephouse=1 if demh==0
replace rephouse=0 if demh==1

*This next variable is the real growth in government agencies
generate realnum=counleg+counpma-terminat

*Now, having selected the lognormal based upon the AIC value, I will continue
*estimating these models.

*Now the parametric models.
*streg, dist (gamma)
*test e(kappa)=1
*test e(sigma)=1
*predict cs, csnell
*stset cs, failure(estate)
*sts generate km=s
*generate H=-ln(km)

*graph H cs cs, c(ll) s(..) t1(Gamma) l1(-ln(Kaplan-Meier)) saving(gamma)
*more

*stset cumyear estate, id(agencyid) 


*streg, dist (gompertz)
*predict csg, csnell
*stset csg, failure(estate)
*sts generate kmg=s
*generate Hg=-ln(kmg)
*graph Hg csg csg, c(ll) s(..) t1(Gompertz) l1(-ln(Kaplan-Meier)) saving (gompertz)

*stset cumyear estate, id(agencyid) 

*streg, dist (loglogistic)
*predict csl, csnell
*stset csl, failure(estate)
*sts generate kml=s
*generate Hl=-ln(kml)
*graph Hl csl csl, c(ll) s(..) t1(Log Logistic) l1(-ln(Kaplan-Meier)) saving(loglog)

*stset cumyear estate, id(agencyid)

*streg, dist (lognormal)
*predict cslg, csnell
*stset cslg, failure(estate)
*sts generate kmlg=s
*generate Hlg=-ln(kmlg)
*graph Hlg cslg cslg, c(ll) s(..) t1(Log Normal) l1(-ln(Kaplan-Meier)) saving (lognorm)

*stset cumyear estate, id(agencyid)

*streg, dist (weibull) time
*predict csw, csnell
*stset csw, failure(estate)
*sts generate kmw=s
*generate Hw=-ln(kmw)
*graph Hw csw csw, c(ll) s(..) t1(Weibull) l1(-ln(Kaplan-Meier)) saving (weibull)

*stset cumyear estate, id(agencyid)

*streg, dist (exponential) time
*predict cse, csnell
*stset cse, failure(estate)
*sts generate kme=s
*generate He=-ln(kme)
*graph He cse cse, c(ll) s(..) t1(Exponential) l1(-ln(Kaplan-Meier)) saving (exponential)

*graph using gamma gompertz loglog lognorm, margin (10) t1(Appendix B.  C-S Residuals of Different Models)


*Now I move on DATA ANALYSIS.  I look first at graphs of the survivor function and then move to full models.
*I will first reset the data.
stset cumyear estate, id(agencyid)

*First I will graph the survivor function by STATUTORY creation.  I first need to create some labels.
label define leglabel 0 "Executive Action" 1 "Statute"
label values leg leglabel

sts graph, by(leg) noborder b2(Age of Agency in Years) t1(Created by Statute)
sts test leg
more
*Now the survivor function by exclusion from REORGANIZATION authority.

*Now recognizing differences across time.
*generate reorchad=1 if reorgaf==1|year>1984
*replace reorchad=0 if reorgaf~=1&year<1985

sts graph, by(reorchad) noborder b2(Age of Agency in Years) t1(Reorganization Authority)
sts test reorchad
more

*Now the survivor function by BUDGET REVIEW.
sts graph, by(afombbud) noborder b2(Age of Agency in Years) t1(OMB Review)
sts test afombbud
more

*Now I will graph the survivor function by the presence or absence of FIXED TERMS.
sts graph, by(term) noborder b2(" ") t1(Fixed Terms)
sts test term
more

*Now I will graph the survivor function by the presence or absence of 
*PARTY BALANCING limitations on appointments.
sts graph, by(pbalanci) noborder b2(Age of Agency in Years) t1(Party Balancing)
sts test pbalanci
more

*After establishing the correct model, I now turn to estimating it.

generate outreorgaf=1 if inreorgaf==0
replace outreorgaf=0 if inreorgaf==1

generate outombrevaf=1 if inombrevaf==0
replace outombrevaf=0 if inombrevaf==1

streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) 
more
*mfx compute
stcurve, survival at1(leg=1) at2(leg=0) xlabel(0(5)50) ylabel(0(.10)1) t1("Statute") b2(Age of Agency in Years) c(ll) s(..)
more
stcurve, survival at1(outreorgaf=1) at2(outreorgaf=0) xlabel(0(5)50) ylabel(0(.10)1) t1("Reorganization Authority ") b2(Age of Agency in Years)c(ll) s(..)

more
stcurve, survival at1(outombrevaf=1) at2(outombrevaf=0) xlabel(0(5)50) ylabel(0(.10)1) t1("OMB Review") b2(Age of Agency in Years) c(ll) s(..)

more
stcurve, survival at1(fixedterm=1) at2(fixedterm=0) xlabel(0(5)50) ylabel(0(.10)1) t1("Fixed Terms") b2(Age of Agency in Years) c(ll) s(..)

more
stcurve, survival at1(partybalance=1) at2(partybalance=0) xlabel(0(5)50) ylabel(0(.10)1) t1("Party Balancing Limitations") b2(Age of Agency in Years) c(ll) s(..)

more
stcurve, survival at1(warrenreg=1) at2(warrenreg=0) xlabel(0(5)50) ylabel(0(.10)1) t1("Regulatory Agency") b2(Age of Agency in Years) c(ll) s(..)

*Now with different measure of preferences.

streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial abdistnc com line temp, dist (lognormal) 
*mfx compute
more
*Now with measures of majority uncertainty.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided bmajorh com line temp, dist (lognormal) 
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided blength com line temp, dist (lognormal) 
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided bseatren com line temp, dist (lognormal) 
more
*Now with cohort indicators.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp y1946 y1947 y1948 y1949 y1950 y1951 y1952 y1953 y1954 y1955 y1956 y1957 y1958 y1959 y1960 y1961 y1962 y1963 y1964 y1965 y1966 y1967 y1968 y1969 y1970 y1971 y1972 y1973 y1974 y1975 y1976 y1977 y1978 y1979 y1980 y1981 y1982 y1983 y1984 y1985 y1986 y1987 y1988 y1989 y1990 y1991 y1992 y1993 y1994 y1995 y1996, dist (lognormal) 
more
*Now with the log of budget size rather than line.
generate lnbudg92=ln(budget)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com lnbudg92 temp, dist (lognormal) 
more
*Now without temporary agencies.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line if temp==0, dist (lognormal) 
more
*Now without the robust estimates of variance.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) robust

streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) cluster (agencyid)

more
*Now with cox model

stcox leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp 
stcox leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp unfrpres unfrmaj unified partintr

*Now I drop two other court-type agencies and see if it affects the results.
drop if agencyid==10458
drop if agencyid==10152
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) 

*Now with the factor score.
factor leg outombrevaf fixedterm partybalance 
rotate
score f1

streg f1 outreorgaf warrenreg corp judicial divided com line temp, dist (lognormal) 

*Now checking for heterogeneity.  First unobserved gamma heterogeneity
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) frailty(gamma)
more
*Now assume some differences in agencies we know about affect the shape of the hazard rate.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) ancillary(temp)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) ancillary(leg)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) ancillary(fixedterm)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) ancillary(partybalance)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) ancillary(outombrevaf)
more
*Now I duplicate this with less assumptions about the hazard.  It could be a totally different baseline hazard
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(temp)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(leg)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(fixedterm)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(partybalance)
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(outombrevaf)


*Now followed with a graph.
streg leg outreorgaf outombrevaf fixedterm partybalance warrenreg corp judicial divided com line temp, dist (lognormal) strata(partybalance)

replace leg=1
replace outreorgaf=0
replace outombrevaf=0
replace fixedterm=0
replace partybalance=1
replace warrenreg=0
replace corp=0
replace judicial=0
replace divided=1
replace com=0
replace line=1
replace temp=0
predict h0, csurv
replace partybalance=0
predict h1, csurv
label var h0 "Party Balancing"
label var h1 "No Party Balancing"
graph h0 h1 _t, c(ll) sort l1title ("Survivor Function by Party Balancing")

